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Abstract 

We present the preliminary results from the search of the non-conservation of lepton flavor number 
in the decay of a r to a lighter mass lepton and a pseudo-scalar meson, performed using e + e~ — > 
t + t~ events collected at a center-of-mass energy near 10.58 GeV with the BABAR detector at the 
SLAC PEP-II e + e~ storage ring. No evidence of such a signal has been found in the data sample 
corresponding to a luminosity of 314.5 fb _1 , and we set an upper limit of 1.6x 10~ 7 at 90% confidence 
level on the decay of r~ — ► fj,~rj. 
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1 INTRODUCTION 



The search for lepton flavor violation (LFV) in charged lepton processes is one of the cleanest ways 
to look for new physics beyond the Standard Model (SM). In the SM with the massless neutrinos, 
the total lepton number and lepton flavors are conserved, so that processes like r~ — > \x~r\ [Q are 
strictly forbidden. The recent discovery of neutrino oscillations 2 indicates that the neutrinos 
have a very small and non-zero mass, and the lepton flavor is not conserved in the neutral lepton 
sector. However, in the simplest extension of the SM, which accommodates the neutrino mass and 
mixing, the expected branching ratios for lepton flavor violating processes in the charged lepton 
sector are too small to be experimentally accessible. Thus, any observation of LFV in the charged 
lepton sector would be an unambiguous signature of new physics beyond the SM [3j. The decay 
r~ — > fj,~r) is of particular interest because it could be enhanced in supersymmetric models [I] 
due to the potentially large coupling of the Higgs boson to the ss pairs and its associated color 
factors. The most stringent upper limit published to date on this decay is Bjjl < 1-5 x at 90% 
confidence level (c.l.) with 154 fb -1 of e + e~ annihilation data collected by the Belle experiment [3]. 

2 THE BABAR DETECTOR AND DATASET 

The results presented in this paper are based upon data collected by the BABAR detector at the 
PEP-II storage ring. Details of the detector are described elsewhere j^]. Charged particles are 
reconstructed as tracks with a 5-layer silicon vertex tracker (SVT) and a 40-layer drift chamber 
(DCH) inside a 1.5-T solenoidal magnet. An electromagnetic calorimeter (EMC) consisting of 6580 
CsI(Tl) crystals is used to identify electrons and photons. A ring-imaging Cherenkov detector 
(DIRC) is used to identify charged pions and kaons and provides additional electron identification 
information. The flux return of the solenoid, instrumented with resistive plate chambers (IFR) and 
limited streamer tubes (LST), is used to identify muons. 

The data sample consists of an integrated luminosity of L= 314.5 fb -1 collected at a center-of- 
mass (CM.) energy y/s near 10.58 GeV. With an average cross section of a e + e -_> T + T - = (0.89±0.02) 
nb [7] as determined using the KK2F Monte Carlo (MC) generator jS], this corresponds to a data 
sample of 279.8 x 10 6 r-pair events. 

The signal MC is generated with a two-body decay r~ — > fj,~r) model, incorporated into the KK2F 
generator employing the TAUOLA decay package [5]. The potential backgrounds including e + e~ — > 
fi + ,t + t~ , (uu,dd, ss mixture), and cc processes were studied using the KK2F, EVTGEN JOJ & 
JETSET generators. We also study bb Monte Carlo events but find that no events pass our 
selection criteria. The detector response is modeled using the GEANT4 simulation package |T2] . 

3 ANALYSIS METHOD 

3.1 Reconstruction and event selection 

We reconstruct the signal process e + e~ — ► t + t~ with one r~ decaying to /i - ??, where r\ — ► 77 or 
7T + 7r _ 7r , using events with zero total charge and with two or four well-reconstructed tracks, where 
none of the tracks originate from conversion of photons in the material of the detector. The event 
is divided into hemispheres by the plane perpendicular to the thrust axis [Ej, which characterizes 
the direction of maximum energy flow in the event, and is calculated using the energy-momentum 
information from all the tracks observed in the DCH and energy deposits in the EMC. 
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For the decay r~ — ► [i~T] (77 — > 77), the signal-side hemisphere is required to contain a pair 
of energy deposits in the EMC each consistent with being a photon above a 100 MeV energy 
threshold |14j and only one track (1-prong hemisphere). Extra energy deposits in the EMC greater 
than the 100 MeV threshold unassociated with any of the identified r daughters are not allowed in 
the signal-side hemisphere. 

For the 3-prong decay r~ — > n~r) (77 — ► 7r + 7r _ 7r°), a pair of energy deposits in the EMC each 
consistent with being a photon above a 50 MeV energy threshold is used to form the tt° candidate. 
We then form two possible 77 candidates by combining the 7r° and the positively charged track 
with either of the two negatively charged tracks in the signal-side hemisphere. The 77 candidate 
closest to its nominal mass (547.75±0.12MeV/c 2 [IE]) is assigned to be the true candidate, and the 
negatively track not associated with the 77, is assigned to be the signal muon candidate. 

We require all possible 77 candidates, constructed out of their respective daughters, to have an 
energy greater than 1.4 GeV and 1.2 GeV for the decay 77 — > 77 and 77 — > it + it~it , respectively. 
We then form the complete r decay chain, and finally accept an event in either one of the decay 
modes of 77 exclusively, based upon the closeness of the signal r candidate to its nominal mass 
(1.777 GeV/ c 2 QH). 

For the 1-prong (77 — * 77) signal candidates, we assign the origin of the two photons from the 77 
to come from the point of closest approach of the signal lepton track to the e + e~ collision axis. For 
the 3-prong (77 — > tt^tt^tt ) decay, the origin of the photons from the 7T° is taken from the common 
vertex found by fitting the three tracks in the signal hemisphere. 

The respective mass windows for the reconstructed 77 candidates are: 

• 0.515 GeV/ c 2 < M(t7 -> 77) < 0.565 GeV/ c 2 , 

• 0.115 GeV/ c 2 < M(tt° 77) < 0.150 GeV/ c 2 , 

• 0.537 GeV/ c 2 < M(t7 -> 71+71-71°) < 0.558 GeV/ c 2 . 

Finally, we kinematically fit for the momentum of the 77 and the 7r° daughters after applying 
77 and 7T° mass constraints on their respective daughters. We then combine the signal muon track 
and the 77 candidates to form t~ candidates. The signal r~ — * 7^77 decays are identified by two 
kinematic variables: the beam-energy constrained r mass (ttt-ec) and AE = E^ + E^ — y/s/2, where 
E^ and E^ are the energy of the \i and the 77 candidate in the CM. frame. These two variables are 
independent apart from small correlations arising from initial and final state radiation. 

The mean and standard deviation of the 777.ec an d AE distributions near their statistical maxima 
for the reconstructed r~ — ► [i~7] (rj — ► 77) MC signal events are: (ttt-ec) = 1777.9 MeV/ c 2 , <r(mEc) 
= 8.25 MeV/ c 2 , (AE) = -13.4 MeV, a{AE) = 40.8 MeV, whereas for t~ -► yrij (77 -> vr+vr-Tr ) 
MC signal events the corresponding quantities are: (tttec) = 1777.7 MeV/ c 2 , 0"(t7Tec) = 5.6 MeV/ c 2 , 
(AE) = —7.1 MeV, a(AE) = 31.0 MeV. The shift from zero in (AE) comes from the leakage of 
the measured photon energy from deposits in the EMC. We do not look at the data events within 
a ± 3d rectangular region centered at ((tttec)> (AE)) in the tttec vs - AE plane, until completing 
all optimization and systematic studies of the selection criteria. 

The r candidate is considered for further analysis if the signal lepton track is consistent with 
being identified as a muon, and not an electron or a kaon, using a set of particle identification (PID) 
criteria using DCH, EMC, DIRC and IFR information, which has a muon identification efficiency 
of about 85% and 65% inside the polar angle coverages of [17°, 57°] and [57°, 155°], respectively, 
for tracks with momentum greater than 0.5GeV/c. This provides the correct association of the 
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reconstructed particles having the r as their mother for 99.8% of selected signal MC events inside 
a ± 5a rectangular region in toec vs. AE plane for both the decay channels. 

To reduce backgrounds from e + e~ — ► p + p~ processes, we require the total CM. momentum of 
all the tracks observed in the DCH and unassociated energy deposits in the EMC on the tag-side 
hemisphere (defined to be the one opposite to the signal-side hemisphere) to be less than 4.75 GeV/c 
for both the channels. 

Other non-T backgrounds are suppressed by requiring the polar angle of the missing momentum 
(Omiss) associated with the neutrinos in the event to lie within the detector acceptance (—0.76 < 
cos Omiss < 0.92), and the scaled missing CM. transverse momentum relative to the beam axis 
(Pmiss/V^) to be greater than 0.082. 

A tag-side hemisphere containing a single track is classified as e— tag, p— tag or h— tag if the 
track is exclusively identified as an electron (e— tag), as a muon (/j— tag), or as neither (h— tag) and 
if the total neutral CM. energy in the hemisphere is no more than 200 MeV. If the total neutral 
CM. energy in the hemisphere is more than 200 MeV and the track is neither an electron nor a 
muon it is classified as p— tag. If the tag-side contains three tracks, the event is classified as a 
3h-tag for the r~ — ► p~r\ (rj — > 77) search, whereas for the 3-prong signal decay channel of 77 we 
use only 1-prong decays of the r in the tag-side hemisphere. 

We assume that the r~ — * p~rj decay fully reconstructs the direction of the r + in the hemisphere 
opposite the signal candidate. We then calculate the invariant mass squared (m 2 ) of the unobserved 
particle, assumed to be v(s), from the r + decay, assuming that the energy of the t + is given by 
the beam energy. 

A cut of — 1.0GeV 2 /c 4 < m 2 < 2.5GeV 2 /c 4 is applied to the e— tag, p— tag and h— tag in the 
search for the r~ — > p~~rj (rj — > 77) decay. For the p— tag, we require — 1.0GeV 2 /c 4 < m 2 < 
1.0GeV 2 /c 4 and for the 3h— tag, we require — 0.2GeV 2 /c 4 < m 2 < 1.0GeV 2 /c 4 in the 77 — > 77 decay 
mode. For 77 — > 7r + 7r~7r° mode, we require — 1.25GeV 2 /c 4 < m 2 < 2.5GeV 2 /c 4 for all tags. In 
addition to the above requirements, we apply cuts on the invariant mass on the tag-side (M tag ), 
calculated from the tracks observed in the DCH and unassociated neutral energy deposits in the 
EMC, to be less than 0.4 GeV/c 2 for the h-tag and 0.6 GeV/c 2 < M tag < 1.35 GeV/c 2 for the p-tag. 

The cuts on the different variables have been applied so as to minimize the expected upper 
limit at 90% c.l. 17 on the branching ratio of the signal r~ —* p~rj decay, in a background only 
hypothesis, estimated using Monte Carlo (MC) simulation of background events. 

After this selection, 10.08% and 5.43% of the r~ — ► pTr\ (77 — ► 77) and r~ — > p~rj (77 — > 7r + 7r~7r°) 
MC signal events survive within a Grand Signal Box (GSB) region defined as: mEC £ [1.5,2.0] 
GeV/c 2 , AE G [-0.8,0.4] GeV. The data distribution of m EC and AE inside the GSB is plotted 
as dots in Figure H About 70% of the selected signal MC events inside the GSB lie within a ±2cr 
rectangular region in the m-Ec vs - plane. These MC events are shown by the shaded region in 
Figure [IJ 

3.2 Background estimation 

The GSB regions excluding the ±3cr blinded region contain 65 and 23 data events, while the 
luminosity-normalized sum of the MC backgrounds yield (74.6 ± 9.0) and (34.3 ± 5.2) events for 
the 2 channels respectively. For the r~ — > p~rj (rj — > 77) channel, 24%, 62% and 9% of the 
MC background events have a true p only, a true 77 candidate only and both true p and true 77 
candidates, respectively. For the r~ — ► p~r] (rj — > 7r + 7r _ 7r°) channel, 4%, 37% and 4% of the 
MC background events have a true p only, a true 77 candidate only and both true p and true 77 
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Figure 1: The AE and toec ditributions for the data (shown as dots) after all the selection criteria 
have been applied, in the search for the — > pTf) (r/ — > 77) decay (left) and the t _ — > /x~r/ (7/ — > 
7T + 7r _ 7r°) decay (right). The shaded region contains 70% of the selected signal MC events inside 
the Grand Signal Box for both modes. The ±2cr signal box is also overlayed. 



candidates, respectively. The major source of backgrounds come from either a mis-identification of 
a pion track as a muon candidate or the in-efficiency in r\ and tt° reconstruction. 

The number of expected backgrounds in the signal box is extracted from an un-binned maximum 
likelihood fit to the distributions of the ttiec an d AE variables for the data events inside the non- 
blinded parts of the GSB region, while the shape of the distributions have been obtained from MC. 
The number of background events {N^ 1 " 1 ) inside the ±2<r signal box is estimated as: 

■\jdata _ llu PDFtot AT data 

JV 2cr — 7 E>nz? X ^GSB-3a 

JGSB-3* ^ Ut tot 

where J 2a PDF to t and /css-3o- PDF to t are the probability density functions (PDFs) integrated 
over the the signal box and the non-blinded parts of the GSB regions, -^gsI?-3o- are t ne number of 
data events in the non-blinded parts of the GSB region, and the PDFtot is defined as: 

PDFtot = (U x PDF ») + Ur x PDF T ) + ([1 - /„ - f T ] x PDF uds ) 

where the and f T are the fractions of and t + t~ background contributions. The PDF^, 

PDF T and PDF U( i s are non-parametric PDFs ^H] obtained from the respective background MC 
events after applying the same final selection criteria. The corresponding projections of the PDFs 
obtained from MC for the individual components and the total fit to the data in the GSB region 
excluding the ±3cr blinded region for mEC an d AE variables are shown in Figure |2 along with the 
data points. The ±2cr signal box is indicated by hatches. 

The number of background events from the fit to the data in the non-blinded parts of the 
GSB region for r _ — > yTr] (r\ — > 77) and r~ — > [i~7] (77 — > 7r + 7r _ 7r°) searches are (0.64±0.08) 
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and (0.07±0.02) events, respectively. The observed and the expected number of background events 
inside the neighbouring boxes in the non-blinded parts of the GSB are shown in Tableland Tabled 
for the respective channels. 

As a cross-check, we also estimate the background assuming a uniform distribution over the 
GSB in rriEC vs - AE 1 plane. This simple interpolation gives the number of backgrounds in signal 
box to be (0.60±0.07) and (0.10±0.02) for t~ -> pTr} (77 -► 77) and r~ -► pTr) (r/ -> vt+tt-tt ), 
respectively. 



r -»■ p t) (?? -► 77) 


(5 - 3)a 


(7-5)a 


(9 - 7)a 


(11 - 9)a 


(11 - 3)<7 


# of observed events 

# of expected events 


4 

2.5±0.3 


3 

3.7±0.5 


5 

4.8±0.6 


7 

5.7±0.7 


19 

16.7±2.1 



Table 1: The number of observed and expected data events inside the 5a — 3a, 7a — 5a, 9a — 7a, 
11a — 9a and llcr — 3d neighboring boxes in the toec vs - AE 1 for the search for the decay — > 
(»7 -» 77)- 



r~ — > (77 — > 7r + 7r~7r°) 


(5 - 3)cr 


(7-5)a 


(9 - 7)a 


(11 - 9)a 


(11 - 3)fJ 


# of observed events 

# of expected events 


1 

0.3±0.1 




0.4±0.1 




0.5±0.1 


2 

0.7±0.1 


3 

1.9±0.4 



Table 2: The number of observed and expected data events inside the 5a — 3a, 7a — 5a, 9a — 7a, 
11a — 9a and 11a — 3a neighboring boxes in the m-Ec vs - A.E plane for the search for the decay 
— > /j,~r) (rj — > 7r + 7r _ 7r°). 



4 SYSTEMATIC STUDIES 

Table 01 is a summary of the relative systematic uncertainties in the signal efficiency estimation. 
The largest systematic uncertainities are due to the signal track momentum and the photon energy 
scale and resolution, which are the tracking and calorimetry errors introduced by applying the final 
cut on the ±2cr signal region. In order to assess the systematic errors associated with these scale 
and resolution uncertainties, the peak and resolution of toec as well as the AE are varied. The 
errors associated with the modeling of each election variable is estimated from the relative change 
in signal efficiency when varying the cut by the Data-MC difference in the mean of that variable. 

Other sources of systematic uncertainities include those arising from trigger and filter efficien- 
cies, tracking and neutral energy reconstruction efficiencies, the PID error associated with signal 
muon track identification, beam energy scale and spread, luminosity estimation and cross-section 
of the e + e~ — > t + t~ process. 

As we use 2.3 million MC signal events, the contribution to the uncertainty arising from signal 
MC statistics is negligible. 

All contributions to the systematic uncertainties are added in quadrature to give a total relative 
systematic uncertainty of 7.5% and 8.6% for the r~ — > fj>~rj (rj — ► 77) and the r~ — ► pT~r\ (77 — ► 
7T + 7T~7r ) channels respectively. 
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Figure 2: The fit to the data distributions of toec an d AE variables in the GSB regions excluding 
the ± 3a blinded region in the search for the r~ — > [i~r] (r/ — > 77) decay (top row) and for the 
r~ — > fi~r] (7/ — > 7T + 7r _ 7r ) decay (bottom row) are shown along with the data distributions (as dots) 
and the background MC components (/i + /U~ in dashed, t + t~ in dotted and uds in dashed-dotted). 
The shaded region is the ±2<r signal box for both modes. 



The dominant contribution to the uncertainity of background estimation arises from the statis- 
tical error on the number of data events surviving the final selection inside the non-blinded parts 
of the GSB, and from the variation of the fitted /„ and f T within ±ler from the fit to the data. 
A small contribution of 1.9% and 0.9% on the relative uncertainities from a bias in the fit arising 
from blinding of the ±3<r rectangular region in the msc vs - plane during estimation of PDFs 
from the MC background samples has been included in the above-mentioned estimate. 
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3 D 


2 7 


Modeling or the selection variables 


3.9 


3.4 


Signal track and photon energy scale and resolution 


3.7 


5.8 


Beam energy scale and spread 


0.4 


0.5 


Luminosity and t + t~ production cross-section 


2.3 


2.3 


Total 


7.5 


8.6 



Table 3: Summary of the relative systematic uncertainties (in %) of the signal efficiency for both 
decay modes t~ — ► \x~~rj {rj — ► 77) and r~ — ► /x~r/ (77 — > 7r + 7r~7r ). 

5 RESULTS 

The upper limit of r~ — > 77^77 is calculated using = N^ L /{2La TT Be), where iV^ is the 90% 
c.l. upper limit on the number of signal events expected within the ±2a box in mEC vs - 
plane, e is the reconstruction efficiency and B is the branching ratio of the 77 decay modes under 
consideration. 

To obtain a combined upper limit, we add the signal efficiencies for the individual channels 
weighted by their respective branching ratios using the formula: 

Be = (Bi x e\ + B 2 x e 2 ) 

where E\ = (7.03±0.53)%, £2 = (3.67±0.32)% are the signal reconstruction efficiencies for the two 
decay modes r~ — > pTr] (77 — ► 77) and r~ — ► ji^rj {rj — ► 7r + 7r~7r ), and 25*i, £>2 are the branching 
ratios of rj — » 77 and 77 — ► 7r + 7r _ 7r (7r — ► 77), respectively. Thus, for the decay modes combined, 
taking into account the errors from similar sources as correlated, we have a combined efficiency for 
reconstructing r~ — > 77^77 of Be = (3.59 ± 0.41)%. This corresponds to a 11.4% relative systematic 
uncertainty on the total selection efficiency. 

We sum the backgrounds inside the ±2<r signal box from the above two channels to obtain a 
background expectation of (0.71 ± 0.08) events for the combined upper limit calculation. Inside 
the ±2<7 signal box, we observe only one event in the search for the decay r _ — * ji~rj {rj — * 77), 
and none in the search of the decay — > ji^rj {rj — ► 7r + 7r~7r°). 

The limit is calculated including all uncertainties using the technique of Cousins and High- 
land 19 j following the implementation of Barlow [201 - I n this technique, MC samples are generated 
according to a Poisson distribution with mean (s + b) where the background, b, and signal, s, are 
each drawn randomly from Gaussian distributions describing their respective PDFs. The values 
of the mean and standard deviation for the background Gaussian are 0.71 and 0.08 events, re- 
spectively. The mean of the signal Gaussian is {2La TT BuLe) an d the standard deviation is the 
error on {2La TT BuLe). The branching ratio, Bjjl, is varied until we find a value for which 10% of 
the sample yields a number of events less than the one event observed in the data. At 90% c.l. 
this procedure gives an observed upper limit for B{t~ — > n~rj) to be 1.6 xlO -7 including effects 
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of systematic uncertainties on the signal efficiency and background estimation for both the decay 
modes combined. Averaging the number of observed events along with its poisson fluctuation, the 
expected upper limit for £>(t~ — > fi~rj) at 90% c.l. is 1.4xl0 - '. 

The r] branching ratio, the selection efficiency, the number of expected and observed background 
events and the observed upper limit at 90% c.l. with the systematic uncertainities are shown in 
Table |U for the two modes separately and combined. 



Decay modes 


Branching Ratio 

of 7?(%) 


Efficiency 

(%) 


Background events 


Upper Limit 
(@90% c.l.) 


Expected 


Observed 


t~ -» iTf] (77 -> 77) 


39.42±0.26 


7.03±0.53 


0.64±0.08 


1 


2.1xl0~ 7 


t~ — > [i~r\ (77 — ► 7r + 7r _ 7r°) 


22.32±0.40 


3.67±0.32 


0.07±0.02 





4.9xl0~ 7 




3.59±0.^ 


11 


0.71±0.08 


1 


1.6xl0~ 7 



Table 4: The 77 branching ratio, the selection efficiency, the number of expected and observed 
background events and the observed upper limit for the two modes separately and combined with 
systematic uncertainity. 



6 SUMMARY 

The search for the SM forbidden decay r~ — > \i~ 77 is performed using 314.5 ftT 1 data. We find 
one event consistent with the signal signature for an expected background of (0.71±0.08) events. 
A 90% c.l. upper limit on the branching ratio is calculated to be 1.6xl0~ 7 including effects of 
systematic uncertainties on the signal efficiency and background estimation. 

7 ACKNOWLEDGMENTS 

We are grateful for the extraordinary contributions of our PEP-II colleagues in achieving the 
excellent luminosity and machine conditions that have made this work possible. The success of 
this project also relies critically on the expertise and dedication of the computing organizations 
that support BABAR. The collaborating institutions wish to thank SLAC for its support and the 
kind hospitality extended to them. This work is supported by the US Department of Energy and 
National Science Foundation, the Natural Sciences and Engineering Research Council (Canada), 
Institute of High Energy Physics (China), the Commissariat a l'Energie Atomique and Institut 
National de Physique Nucleaire et de Physique des Particules (France), the Bundesministerium fur 
Bildung und Forschung and Deutsche Forschungsgemeinschaft (Germany), the Istituto Nazionale di 
Fisica Nucleare (Italy), the Foundation for Fundamental Research on Matter (The Netherlands), the 
Research Council of Norway, the Ministry of Science and Technology of the Russian Federation, and 
the Particle Physics and Astronomy Research Council (United Kingdom) . Individuals have received 
support from the Marie-Curie IEF program (European Union) and the A. P. Sloan Foundation. 



15 



References 

[1] Charge conjugate mode for all the decays are implied throughout this paper unless otherwise 
stated. 

[2] B. T. Cleveland et al, Astrophys. J. 496, 505 (1998); Y. Fukuda et al. [Super-Kamiokande 
Collab.], Phys. Rev. Lett. 81, 1562 (1998); Q. R. Ahmad et al. [SNO Collab.], Phys. Rev. 
Lett. 89, 011301 (2002); M. H. Ahn et al. [K2K Collab.], Phys. Rev. Lett. 90, 041801 (2003); 
K. Eguchi et al. [KamLAND Collab.], Phys. Rev. Lett. 90, 021802 (2003). 

[3] E. Ma, Nucl. Phys. Proc. Suppl. 123, 125 (2003). 

[4] M. Sher et al, Phys. Rev. D66, 057301 (2002). 

[5] Y. Enari et al. (Belle Collaboration), Phys. Lett. B622, 218 (2005). 

[6] B. Aubert et al. (BABAR Collaboration), Nucl. Instrum. Methods Phys. Res., Sect. A 479, 1 
(2002). 

[7] The uncertainty was conservatively estimated by comparisons of the cross section between the 
generators KK2F |5] and K0RALB: S. Jadach and Z. Was, Comput. Phys. Commun. 85, 453 
(1995). 

[8] B. F. Ward, S. Jadach, and Z. Was, Nucl. Phys. Proc. Suppl. 116, 73 (2003). 

[9] S. Jadach, Z. Was, R. Decker, and J. H. Kuhn, Comput. Phys. Commun. 76, 361 (1993). 
[10] D. J. Lange, Nucl. Instrum. Methods Phys. Res., Sect. A 462, 152 (2001). 
[11] T. Sjostrand, Comput. Phys. Commun. 82, 74 (1994). 

[12] S. Agostinelli et al. (GEANT4 Collaboration), Nucl. Instrum. Methods Phys. Res., Sect. A 
506, 250 (2003). 

[13] S. Brandt et al, Phys. Lett. 12 57 (1964); E. Fahri, Phys. Rev. Lett. 39 1587 (1977). 

[14] All measured quantities are quoted in the laboratory frame, unless otherwise mentioned. 

[15] S. Eidelman et al. Phys. Lett. B 592, 1 (2004). 

[16] J. Z. Bai et al. (BES Collaboration), Phys. Rev. D 53, 20 (1996). 

[17] G. J. Feldman and R. D. Cousins, Phys. Rev. D 57, 3873 (1998). 

[18] K. Cranmer, Comput. Phys. Commun. 136, 198 (2001). 

[19] R. D. Cousins and V. L. Highland, Nucl. Instrum. Methods Phys. Res., Sect. A 320, 331 
(1992). 

[20] R. Barlow, Comput. Phys. Comm. 149, 97 (2002). 



16 



